Soft modes near the buckling transition of icosahedral shells 
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Icosahedral shells undergo a buckling transition as the ratio of Young's modulus to 
bending stiffness increases. Strong bending stiffness favors smooth, nearly spherical 
shapes, while weak bending stiffness leads to a sharply faceted icosahedral shape. 
Based on the phonon spectrum of a simplified mass-and-spring model of the shell, 
we interpret the transition from smooth to faceted as a soft-mode transition. In 
contrast to the case of a disclinated planar network where the transition is sharply 
defined, the mean curvature of the sphere smooths the transitition. We define elastic 
susceptibilities as the response to forces applied at vertices, edges and faces of an 
icosahedron. At the soft-mode transition the vertex susceptibility is the largest, but 
as the shell becomes more faceted the edge and face susceptibilities greatly exceed 
the vertex susceptibility. Limiting behaviors of the susceptibilities are analyzed 
and related to the ridge-scaling behavior of elastic sheets. Our results apply to 
virus capsids, liposomes with crystalline order and other shell-like structures with 



2 



icosahedral symmetry. 
PACS numbers: 



I. INTRODUCTION 



Virus capsids p| and other structures such as coUoidosomes [2[ and hposomes 



m 



consist of thin shells of spherical topology that frequently exhibit icosahedral symmetry. A 



popular simplified model 

mm 

replaces the shell with a triangulated network of masses 
and springs (see Fig. [1]). This network consists of five- and six-coordinated vertices, with 
the five-coordinated vertices aligned with the five-fold icosahedral symmetry axes. Five- 
coordinated vertices may be considered as +27r/6 disclinations within an otherwise six- 



coordinated lattice. 

spherical shells ay 



'hese disclinations are absent in conventional continuum models of 



Elastic properties of the capsid can be mimicked by suitably adjusting the spring constants 
to obtain the desired Young's modulus Y and by imposing a curvature energy to obtain 
the bending modulus k. Strains associated with the disclinations cause the network to 
buckle 12|, transforming the shape from smooth and nearly spherical to strongly faceted 




FIG. 1: Triangulated network oi P = Q = 2. Colors identify local environments with 5- fold vertices 
shown in blue. 



FIG. 2: Shell shape above the buckling transition for P = 128, Q = shell with kg = 1 and kh = 16 
yielding Foppl-von Karman number 7=930. Color coding is logarithmic according to total elastic 
energy (violet=low, red=high) 



and nearly icosahedral [5]. A dimensionless parameter controls the transformation. We 
define the Foppl-von Karman number 



7 



K 



(1) 



where i? is a linear dimension of the shell. The buckling occurs when 7 exceeds a valu e 7^, of 
order 10^ (see Fig. [2]). For the virus HK97, which appears to facet as it matures 13|, [ij], 7 



reaches a value of order lO'' according to the estimate of Ref. Q|. Varying the pH of solution 



can alter 7, with the range 100-900 reported for the virus CCMV 15|, ll6|]. At much larger 



crysta 



20 



21 



line order, an 



22|. 



values of 7 (in excess of 10^) which should characterize liposomes with^ 
interesting phenomenon known as "ridge scaling" emerges Q, Q, 

Caspar and Klug jl| classify icosahedral structures by a pair of integers {P,Q). A pair 
of five-coordinated vertices is connected by a path consisting of P edges in some given 
direction and Q edges in a direction 60° to the left (e.g. between two blue vertices via a red 
vertex in Fig. [T]). The T-number of the network, T = P"^ + PQ + gives the number of 
vertices as A^^ = lOT + 2. There are always 12 five-coordinated vertices, so the number of 
six-coordinated vertices is 10 (T — 1). Structures with P and Q both nonzero and P ^ Q 
are chiral, such that (P, Q) and {Q, P) are mirror images. Their symmetry group is the 
60-element icosahedral rotation group Y. Structures with either P or Q = 0, or with P = Q 
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are nonchiral. Their symmetry belongs to the 120-element group Ij, = YxZ2, which should 



23| Y'. 



not be confused with the 120-element icosahedral double group 

We exploit the rotational symmetry group to analyze the normal modes of the network 
model by diagonalizing the Hessian matrix of the elastic energy. Eigenvectors represent 
characteristic modes of deformation, which transform according to irreducible representa- 
tions of Y, and the corresponding eigenvalues measure the mechanical stability. Because the 
buckling occurs in a symmetric fashion, the corresponding modes must exhibit full icosa- 
hedral symmetry. Nondegenerate modes transform as the unit representation. Tracking 
these nondegenerate eigenvalues reveals a softening and also a mixing of modes as 7 passes 
through 7b. 



Other studies consider more microscopis elastic network models [24, l25|] that place nodes 
at every atom in the amino acid chains. These studies find that the displacements 
during maturation (i.e. as the virus goes through the buckling transition) can be accurately 
represented using a superposition of only the lowest few nondegenerate modes, consistent 
with our expectations. 

Section [Til of this paper reviews the continuum-elastic theory for deformations of planes 
and spheres, to establish notation and for comparison with our later numerical results. Our 
network model is defined in section [Till and applied to the special cases of disclination- 
free triangular lattices, single disclinations of positive charge, and icosahedral structures 
of spherical topology containing twelve disclinations. Low-lying eigenvalue spectra reveal a 
sharp buckling transition in the case of a single disclination but a broadly smeared transition 
for the icosahedral case. Following Ref. we find that the positive curvature of the sphere 
plays a symmetry-breaking role analagous to an applied magnetic field at a ferromagnetic 
phase transition. 

The final section (lIVp applies forces to selected points on a plane or a shell to probe the 
elastic response of the network as a whole. The resulting displacements define susceptibilities 
which diverge in the case of the single disclination. In the case of the icosahedron, we find 
the effective stiffness (inverse of the susceptibility) drops most rapidly at 7;, for forces applied 
at five-fold symmetry axes, but the stiffness falls off more rapidly for forces applied at two- 
and three-fold symmetry axes for 7 > 7;,. We analyze these susceptibilities in limiting cases 
of small and large 7. 
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II. CONTINUUM-ELASTIC THEORY 



The general elasticity theory of membranes can be expressed in coordinate-free form ll|, 



261]. Let M be a manifold (a two-dimensional smooth surface embedded in three dimensional 
space) assumed to be in mechanical equilibrium. Now impose tangential deformation u(x) 
and normal deformation ({x) corresponding to displacements of points x on the surface. Let 
Qa/B and Cai3 be the metric and curvature tensors respectively of M after distortion. Greek 
indices take values 1 and 2 corresponding to the dimensions of M. Define the strain tensor 

Uaf3 = Uaf3 + (Cap (2) 

where Ua/s = ^{DaUp + DpUa) and indicates covariant differentiation with respect to Xa. 



The trace measures dilation, while 



Sad = Uaf3 - ^ga/sU^I (3) 



measures shear strain. Bending of M is characterized by mean curvature H = ^C^ and 
Gaussian curvature K = det C. 

The free energy density at x contains dilation, shear and bending contributions, 

/(x) = + + (4) 

and can be integrated over M to obtain the total free energy 

^ /(x)v/dit^rf2x. (5) 

The separate contributions are 

/d(x) = i(A + ^)(f/^^)2 (6) 
/,(x) = fiS'^f'Sap 
/fc(x) = ^k{2H - cof + kgK. 



The elastic constants A and are the Lame constants [27|. The 2D area (bulk) modulus 
B = \ + fi, while fi itself is the shear modulus, and the 2D Young's modulus Y = 4/i(A + 
yu)/(A + 2/i). Upon integration over the surface M, the Gaussian curvature term becomes 
constant, and we neglect this term henceforth. Likewise, we set the spontaneous curvature 
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Co = 0, thus assuming the manifold M would be flat in the absence of constraints associated 
with the spherical topology. Effects of Cq 7^ are discussed in Ref. Q]. 
Given the elastic free energy F we obtain the stress tensor 

^"<^> ^ 

whose divergence yields the tangential force 

F" = D^a''^. (8) 

The normal force is given by 

In mechanical equilibrium the stress tensor and normal force vanish. Slightly out of equilib- 
rium, to flrst order in the displacements, special forms of u(x) and known as normal 
modes solve the eigenvalue equation 

(F,iV) = -A(u,C). (10) 

When displaced from equilibrium according to the k^^ normal mode (u^, (k), the free energy 
increases by 

AF=U, y"(|ufep + C|)rfx. (11) 

According to the equipartition theorem the modes fluctuate with thermal energy AF = 
kBT/2 and amplitude /(|ufcp + Qdx = 2kBT / K^. 

Time dependence of the strains depends on the equations of motion. In the overdamped 
case we write 

«^(x) = rF^(x), c(x) = riv(x) (12) 

where we take F proportional to an inverse viscosity as in a Stokes-Einstein relation. In this 

281 carries out a more 



case a normal mode decays in time with a decay rate uj = FA. Ref 
thorough investigation of flat membranes coupled to fluid flow. In the absence of damping 
we write 

pM^(x) = F^(x), pC(x) = iV(x) (13) 



with p the 2D mass density. A normal mode now oscillates in time at frequency lu = ^J~KJp. 
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A. Deformations of a Plane 

An infinite flat elastic sheet in equilibrium has no curvature, so for small perturbations 
the energy decouples into contributions from the in-plane strain u and perpendicular dis- 
placement C- 

/ = iA(«^)2 + + \K{/\Cf (14) 

Here A = DaD°' = is the usual 2D laplacian operator and V the usual gradient. By 
differentiating the energy we obtain the forces 

F = (A + /i) VV ■ u + fiAn (15) 

and 

N = -kA'^C (16) 

Because the in-plane and out-of-plane displacements decouple, we solve them separately. 
The solutions are based on the plane wave function 

V'k(r) = e^*^-"- (17) 

which is an eigenfunction of the Laplacian operator A^/'klr) = — ^^V^k(i")- In-plane normal 
modes are expressed as longitudinal waves 

UL(r) = V^k(r) = ike''^-' (18) 

and transverse waves 

uy(r) = z xul = i{Ky - kyx)e'^"' . (19) 

Note the identities V x = and V ■ uj- = 0, as expected for longitudinal and transverse 
waves. These waves are eigenvectors of the in-plane force eq. ( ITSl) provided their eigenvalues 
obey the longitudinal and transverse dispersion relations, respectively 

Ai = (A + 2/i)A;2 (20) 

and 

At = /ifcl (21) 

Perpendicular out-of-plane waves 

up(r) = zij^ir) (22) 
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obey eq. (fT6|) subject to the perpendicular wave dispersion relation 

Ap = Kk\ (23) 

For future reference we recast the normal modes in plane-polar coordinates (r, 0), replac- 
ing the plane- wave function ip]i{r) with cylindrical Bessel functions 

^Pkm{rA) = Jm{kr)e'"''^. (24) 
The Laplacian operator takes the form 

A - -— (r—^ + — — (25) 

r dr \ dr J dcj)^ 

Like the plane wave function ^/'k(r), waves of type fl24|) are eigenfunctions of the Lapla- 
cian operator, A'?/'fcm(r, 0) = —k'^ipkm{f',4')- Upon defining normal modes Ul,Ut,Up as in 
eqs. ffTSjl . f|T9|l and fl22l) the longitudinal, transverse and perpendicular dispersion relations 
given in eq. ( l20l) . (I2T1) and (l23i) result. These polar forms generalize nicely to conical and 
spherical geometries. 

B. Deformations of a Sphere 

Now we redo the prior calculation of section [TiAl for the case of small perturbations around 
a sphere of equilibrium radius R. In this case the unperturbed manifold has constant mean 
curvature Hq = 1/R. Consequently the free energy acquires a term coupling the in-plane 
and normal strains through the dilation energy. 

fd = l{X + f^){u:^ + 2C/Ry (26) 

In the above, the Laplacian operator takes the form 

1 d f . d\ 1 d"^ 



Notice we include the integration measure -y/det g along with the bending energy /{,, because 
it contributes the term {DaCY- The v^det g factor is not needed in fd or fs because these 
are already second order in the deformation. 
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Taking functional derivatives of F yields the stress tensor, in-plane and normal force 

a"/^ = Xg^^'iul^ + |C) + 2/x(w"^ + ^g^^Q (28) 
F^=iX + /i)Z}^(«^ + |C) + /i(A + j^)u^ 

iv = -(A + /i)(|«^ + -lc) 



where we define 



C = D^D'^DpD^ + Jj^a^". (29) 



The extra ^u^ /B? in eq. (l28l) for comes from commutation of covariant derivatives. The 
final, second derivative, term in ( l29l) comes from integrating by parts the square of the first 



derivative in /fe-^/det g. 

Take the spherical harmonic Yim{d,4>) as the basic deformation, analagous to the plane 
wave e^^ '^ in eq. f|T7|) or the cylindrical wave Jm{kr)e^^'^ in eq. fl2^ . The spherical harmonic is 
an eigenfunction of A with eigenvalue —l{l + l)/B? and an eigenfunction of C with eigenvalue 
/(/ — !)(/ + !)(/ + 2)/R^. By analogy with the procedure for plane waves in flat space, we 
take derivatives as 

ul = RVYi^ ul = RD'^Yirn (30) 

where e is the alternating tensor. We also deflne 

Up = rYim (31) 

These functions are linear combinations of the "Vector Spherical Harmonics" V/m, and 
X^m, which form a complete set of orthogonal functions for expanding vector flelds on the 



surface of a sphere [23|, |29|] . Notice that the transverse mode is proportional to the angular 
momentum operator acting on Yim, thus identifying it with the vector spherical harmonic 
Xim. The longitudinal and perpendicular modes, ul and up, are linear combinations of V^^ 
and yVim- Note that the longitudinal and transverse modes ul and exist only for / > 1, 
while Up exists for / > 0. 

The transverse mode u^ is divergenceless {u'^ = 0) and hence creates no perpendicular 
force and no longitudinal force (the gradient part of F^). In fact, it is an eigenfunction 
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[n ) 







II 




Mpp j 







of the force (128|) . Upon taking into account the commutation of covariant derivatives, we 
find = [fi{l — !)(/ + 2)/R'^]u^ from which we obtain the eigenvalue 

As expected, A^- = for / = 1 because these modes correspond to rigid rotations. 

In contrast to the transverse modes, the longitudinal and perpendicular modes and 
Up are coupled in both the tangential force and perpendicular force A^. In matrix form. 



(33) 



Setting u^, as in eq. (130|) and setting ( as the radial component of up in eq. (I3T1) . the matrix 
elements of M become 

M,, = (A + /i)^ + ^(^^f±^ 

MiP = (A + /i)f 

Mp, = (A + /.)^ 

Mpp = (A + /i)^ + ^ ^'-W)a+2) 

The eigenvalues of this matrix, A-t, are the desired normal mode eigenvalues A. For the 
special case / = 1 the eigenvalues are A_ = and A+ = 6(A+/i)/-R^. The vanishing eigenvalue 
A_ corresponds to rigid translation (for example, the north and south pole displace upwards 
perpendicular to the shell while the the equator displaces upwards tangent to the shell). The 
finite eigenvalue A+ corresponds to an "optical" mode in which polar and equatorial regions 
displace in opposite directions (for example, the north and south poles displace upwards 
while the equator displaces downwards). 

The spherical solution should go smoothly to the flat space solution in polar coordinates 
as the sphere radius R ^ oo. This correspondence can be verified by holding r = R6, 



k = l/R and m fixed, and noting 30 1 



}^ \ = i-'^rJmikr)e'"''^. (35) 

>oo V Z< + i 

In addition, the eigenvalues should approach their proper limits. Clearly At approaches 
its flat space limit fl2T]) . To check A^ p, note that the eigenvalues A-t of the matrix flM|) 
approach (A + 2yu)/(/ + l)/i?^ and k,{1 — !)/(/ + !)(/ + 2)/ R'^ in the limit of large sphere radius 
i?, yielding the flat space limits Eqs. (l20i) and (l23ll . 
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III. MASS AND SPRING MODEL 

We now introduce the discrete mass and spring model for which numerical calculations 
will be performed. This model is also closer to reality for liposomes and colloidosomes, which 
consist respectively, of discrete lipid molecules and colloidal particles, and also for viruses, 
which consist of an aggregation of discrete protein subunits known as capsomers. Let rj be 
the position of mass i and hj be the orientation of plaquette /. A plaquette is a set of three 
masses joined to each other by springs, and we take the normal in the outward direction. 
Following 5| we define 

ns = ^J2^\r.-T,\-a)' (36) 



2 



and 



m = ^J2\^r-nj\^ (37) 

(/J) 

and set the unstretched spring length a = 1. Here {ij) denote pairs of nearest-neighbor 
vertices, and {IJ) denote pairs of adjacent (edge-sharing) plaquettes. In the continuum 
limit the discrete model reproduces the continuum system with elastic constants 

F = -^fc, K = ^h, (38) 



Foppl-von Karman number 



7 = = 39 

K 6kh 



Lame coefficients and bulk modulus 

A = /i = ^fc, B = y^K, (40) 
and 2D mass density (taking the vertex mass m = 1) 

p = 2/V3. (41) 

A. Deformations from fiat-space 



Consider a regular six-coordinated triangulated network of masses and springs. As before 
we start with the plane wave function (fT7|) and take its gradient to obtain the longitudinal 
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sound wave. The dispersion relation is simplest for wavevector k in the y direction (chosen 
to lie midway between two near-neighbor bonds) 

\/3 

KL = 'iK{l-cos{^kya)) (42) 

Taking the cross product with z yields the transverse sound wave with dispersion relation 

/3 

Ar = fc.(l-cos(^A;j,a)) (43) 

Finally, taking the perpendicular displacements as the planewave yields the perpendicular 
modes with dispersion relation 

a/3 

Ap = h{2 - 2 cos {^kya)f (44) 

In the continuum limit ka « 1 these dispersion relations revert to the prior results of 
continuum elastic theory. 



B. Buckling of a Plane into a Cone 

Upon introducing a five-fold +27r/6 disclination into the flat triangulated network dis- 

n 

cussed previously, strain energy accumulates [12] and grows without bound as the radius 
R of the network increases. At some specific "buckling radius" -Rj, it becomes energetically 
favorable to buckle out of plane, trading a reduction in strain energy for a cost in bending 
energy. The trade-off is measured by the Foppl-von Karman number 7. Small 7 favors fiat 
networks, while larger 7 favors buckling into a conical shape. 

In the following we analyze the vibrational spectrum of the network as it passes from flat 
to conical. Rather than vary the radius, we hold R fixed and vary the bending stiffness. 
Large k^ opposes buckling and the network lies fiat, while small kh allows buckling out of 
plane into a cone. For the network of radius R = 8a analyzed below, buckling occurs for 
kb ~ 0.71. As R increases the threshold value of kb grows as R^ so that 7 approaches the 
limiting value 75 ~ 154 

Eigenvectors of the Hessian matrix form basis functions for representations of the symme- 
try group of the structure [sil . Eigenvectors sharing a common eigenvalue form the basis for 
an irreducible representation. Thus the patterns of degeneracy follow the dimensionalities of 
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FIG. 3: Triangulated network of radius 8a and spring constant kg = 1 with a single 5-fold disclina- 
tion at center, (top) Eigenvalue spectrum color coded according to degeneracy. Note the nonde- 
generate Ix mode that goes to zero at the buckling transition, (bottom) energy, cone height and 
susceptibility. 



14 





m 


ICo 


2C5 


2C| 


5(T„ 


Ai 





1 


1 


1 


1 


A2 





1 


1 


1 


-1 


El 


1 


2 


r-i 


— T 





E2 


2 


2 


— r 








TABLE I: Character table of Cs^. C„ denotes conjugacy class of order n. Values of m denote 
in-plane angular momenta, r = (\/5 + l)/2 is the Golden Mean. 

the irreducible representations, as can be seen in Fig. [HH- Likewise the eigenvectors exhibit 
special symmetry properties associated with subgroups of the full symmetry group. 

The symmetry point group of the cone is C^^ in general, corresponding to five-fold rota- 
tions around an axis passing through the five-coordinated vertex, together with reflections 
in vertical planes passing through this axis (see Table [T]). For the specific case of the fiat net- 
work, the group is even higher, -D5/1, adding reflections in the horizontal plane, and two-fold 
rotations around axes lying within the plane. For both groups all irreducible representations 
are either 1- or 2-dimensional, so all nonzero eigenvalues must be nondegenerate or two-fold 
degenerate. Of course, there must be a sixfold degeneracy of zero eigenvalues, corresponding 
to rigid translations and continuous rotations (not belonging to the finite point group) that 
leave the energy invariant. 

For the group D^^, the irreducible representations are based on those of C^^ supplemented 
with an additional label g, u according to whether they are even {g) or odd {u) under reflec- 
tion through the horizontal plane at- The requirement that each irreducible representation 
be either even or odd under at requires that each mode be polarized either fully in-plane or 
fully perpendicular. 

Let Ai be the lowest nonzero eigenvalue. Its eigenvector ei is polarized strictly perpendic- 
ular to the sheet and transforms as the irreducible representation A2u- Its value is nonzero 
at the origin. The energy of mode i varies as Ajof where measures the amplitude of the 
mode. Mechanical equilibrium thus demands that all eigenvalues (other than the six zero 
modes) be strictly positive. In particular it requires Ai > 0. However, if we monitor the 
value of Ai as a function of 7 (Fig. [3^) we find it crosses through zero at 7f,. 
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For small deformations we express the energy as 

i 

Now set Ai = c(7{, — 7). The mechanically stable minimum energy structure is perfectly 
flat (oj = 0) for 7 < 76, but it buckles out of plane for 7 > 7;,, in a shape described by the 
eigenvector ei, with amplitude growing as 1/7 — 7b. Meanwhile the energy drops as (7 — 76)^- 
These effects can be seen in figure [3b. 

For 7 > 7b, figure [3k shows the spectrum of vibrations around the mechanically stable, 
buckled structure. Note that A'^ (the lowest nondegenerate eigenvalue) becomes positive 
again. 



C. Buckling of Spherical Shells 

1. P = 1,Q = icosahedron 

Table [Til presents the character table of the 60-element icosahedral rotational symmetry 
group Y, which has 5 irreducible representations. The conjugacy classes are labeled Cn, 
where n is the order of elements in the class, so that an element of C„ corresponds to 
a rotation by 2TT/n. Recall that the spherical harmonics Yi^ form basis functions for the 
irreducible representations of the continuous rotation group 5*0(3), and therefore they induce 
representations (in general reducible) of Y. For a given angular momentum / and rotation 
angle 6, the character is 

Irreducible representations of Y are labeled in Table [III according to the lowest angular 
momentum / under which they transform. Of particular interest is the representation Fi 
corresponding to angular momentum / = 1. This is the representation under which three- 
dimensional vectors transform. 

The simple icosahedron has 12 vertices, 20 faces and 30 edges. Since we place masses on 
the vertices, our eigenstates are functions defined only at vertex positions. Arbitrary scalar- 
valued functions can be expressed as linear combinations of the basis functions of the "regular 
representation", one of which is concentrated at each icosahedron vertex. The characters 
Xr of the regular representation equal the number of vertices that remain stationary under 
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TABLE II: Character table of Y. Cn denotes conjugacy class of order n. Values of / denote angular 
momenta. R is the "regular representation" and V the "total vibrational representation" discussed 
in subsection IIII C II 

a given symmetry operation. Our vibrational modes are wector-valued functions on the set 
of vertices and thus can be expressed as linear combinations of the product of the regular 
representation R times the representation Fi corresponding to a three dimensional vector. 
We call the resulting product representation the "total vibrational representation" [23|], and 
its characters xv = XRXF-^■ 

Reducible representations can be decomposed into their irreducible components using 
orthogonality properties of character tables. In particular we obtain the decomposition 



= A © 3Fi © F2 © 2G © ZH. 



(47) 



Of the three occurrences of the vector representation Fi we know that two must correspond 
to rigid global translations and rotations. These leave the energy invariant and hence are 
zero frequency modes. The nondegenerate mode transforming as the unit representation A 
must correspond to a "breathing mode" in which all vertices displace equally in the radial 
direction. We find that the remaining modes have specific interpretations in terms of vector 
spherical harmonics, as listed in table lllli 



2. Higher Order Icosahedra 

As the icosahedron is subdivided and the total number of vertices grows, the classification 
of modes into irreducible representations remains similar, but each irreducible representation 
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Formula Irrep g Comments 



0.00000 2 X Fi 6 Translations + rotations 

0.58579 2 — a/3 Ha 5 Mixed contains and 

0.76393 Vb/R"^ -2 F2 3 Radial contains rYsm 

1.00000 1 H 5 Tangent Xs™ 

1.80901 1 + r/2 Ga 4 Tangent Xg^ 

2.76393 ^/5/i?2 A 1 Radial ryoo breathing mode 

3.00000 3 Fi 3 Mixed Vi^ 

3.41421 2 + ^2 Hb 5 Mixed 

3.42705 1 + 3r/2 Gb 4 Tangent 

TABLE III: Vibrational eigenvalues for P = 1,Q = icosahedron with a = kg = 1, kb = 0. A is 
eigenvalue and g is degeneracy. R = Vl + = 0.95106 is the radius of the icosahedron. 

now occurs many times. Fig shows the lowest frequency modes for a P = 8,Q = 
icosahedron with = 642 vertices. To obtain this figure, we set ks = 1, and varied kb. For 
each value of kb we relaxed the structure to mechanical equilibrium using steepest descents, 
evaluated the Hessian matrix by numerical differentiation, then diagonalized the matrix. 
The Foppl-von Karman number is defined as in eq. ( l39i) . where now R is defined as the 
root-mean-square radius (defining R instead as the mean radius [Is!] has little impact below 
or near the buckling transition and results only in a slight rescaling as 7 grows large) and 
takes values in the range 6.6-7.6 for the P = 8,Q = icosahedron. 

Owing to rotation and translation invariance of the total energy, we always have a 6-fold 
degenerate mode of zero eigenvalue. The remaining eigenvalues fall into the classification of 
icosahedral symmetry introduced in Table HTl 

At low 7, when the shape is spherical in the continuum limit of large radius, and the 
energy cost of bending dominates over the energy cost of stretching or shearing, the lowest 
frequency nondegenerate mode is a "breathing" mode, corresponding to a sphere with oscil- 
lating radius. Perturbing the radius by an amount ( (i.e. adding mode up = f^) increases the 
energy by SnBC^ while displacing vertices by (. Identifying the energy with ^NyAbreathing, 
and noting the area modulus B = we find eigenvalue ^breathing = Sny/Skg/N^ which 

fits well to the data in Fig. HI 
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FIG. 4: Lowest frequency modes of P = 8, Q = icosahedron with = 642 vertices, (top) Color 
coded according to degeneracy, (bottom) Nondegenerate modes only. Arrows locate eigenvalue 
^breathing = SvTV^/A't,. Note the buckling mode (red) dips close to zero near the buckling transition. 
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At higher frequencies, where the wavelength of the modes becomes small compared to the 
radius of curvature, we expect that the eigenvalues should revert to their fiat space limits as 
discussed in section 111 B[ The validity of this hypothesis is demonstrated in the dispersion 
relations shown in fig. [5l Here we plot the vibrational frequencies (i.e. the square roots 
of Hessian eigenvalues) as functions of the equivalent wave number, defined as the angular 
momentum index / divided by the radius R. The radii of the circles represent the projections 
of the eigenvectors onto the vector spherical harmonics X;^ (top); and the longitudinal and 
transverse eigenf unctions u.^ and u^- (middle and bottom). The solid lines are the predictions 
of continuum elastic theory for the plane, eqs ( l20ll22l) . 

Soft-mode behavior at the buckling transition is less pronounced than in the case of the 
cone. The crossover from spherical to faceted shape, which occurs gradually for 7 ~ 100 — 
1000, preserves the icosahedral symmetry. As such, the displacements respect icosahedral 
symmetry. If the transition is due to a "soft mode", this mode itself must be invariant 
under operations of the icosahedral symmetry group. That is, it must transform as the 
unit representation and therefore must be nondegenerate. The soft mode is best seen in 
figure Hb, where only the nondegenerate modes are shown. Always the lowest frequency 
nondegenerate mode is an Z = m = "breathing" mode, and as just discussed its frequency 
does not depend significantly on /c^. However, the next occurrence of the unit representation, 
at Z = 6, contains a mode, of type up and labeled Igm; that does indeed soften and mixes 
with the breathing mode in an avoided crossing around 7 ~ 400. Another I = 6 mode, of 
type ut and labeled X2m, is prevented by symmetry from mixing with the up mode. A 
series of other nondegenerate modes Yim (/ = 10,12,20,...) soften at higher 7 values and 
mix with the other soft modes. 

Around 7?, the buckling mode consists predominantly of / = and / = 6 spherical har- 
monics, with a small admixture of / = 10 and higher harmonics. The weight of this mode is 
concentrated in the vicinities of the icosahedron vertices, and it has strong overlap with the 
displacements of vertices under the buckling transition. 

The forbidden crossing of the buckling and breathing mode smears of the buckling tran- 
sition, because Ahuckiing > ^breathing > 0) prcvcuts the eigenvalue of the buckling mode from 
actually crossing zero. This contrasts with the case of the disclinated fiat sheet buckling 
into a cone, where the eigenvalue does indeed cross zero. For the sheet-to-cone transition 
the analogue of the breathing mode is just a zero energy translation, rather than a finte 
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FIG. 5: Vibrational frequencies plotted versus wave number q = l/R, where / is the angular 
momentum index. Data is for P = 8,Q = icosahedron with kg = l,kb = 0.1, 7 = 653. The radii 
of the circles indicate the sizes of the various projections. 
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frequency radial displacement. Also, up-down symmetry of the plane allows the crossing of 
the buckling mode (which is odd) with this translation. On a sphere the symmetry breaking 
between inside and outside the sphere causes the breathing mode to mix with the buckling. 

Owing to the smearing, the value of 7^ is not uniquely defined for the sphere-to- 
icosahedron transition. Reported values range from 130 to 260 based on fitted energy 
models We observe the avoided crossing around 7 ~ 400. 



IV. SUSCEPTIBILITIES 



A. Cones 



The soft-mode transition is a genuine sharp phase transition for the buckling of a discli- 
nated sheet into a cone. We already discussed the order parameter (height) and energy 
variation through the transition, in section UlI B[ Now we consider the susceptibility, namely 
the response of the order parameter to an applied field. In this case we examine the response 
of the buckling height to a point force applied at the disclination. 

Assume the height of the cone (i.e. the vertical displacement of the 5-coordinated particle 
at the center) is given hj h = Piflj, where again Oj is the amplitude and Pi measures the 
projection of the mode i onto the height variable. Then in the presence of an applied force 
we express the energy as 

E = Y,il^^ci■-FP,a,), (48) 

i 

which differs from (l45l) by the work done against the applied force. Minimizing the energy 
yields = FPi/Ai resulting in total height h = FY^^Pf/Ki and susceptibility 

i 

Thus a vanishing eigenvalue, say Ai, passing linearly through zero at 7 = 7^, translates 
immediately into a diverging susceptibility. This divergence is evident in Fig. [3)d. Note that 
the amplitudes differ on the two sides of 7?, because in one case we perturb around a flat 
network while in the other we perturb around a buckled cone. 
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FIG. 6: Stiffness (inverse susceptibility) to forces applied with icosahedral symmetry at vertices, 
edges and faces. Dashed blue lines show 15/ log (0.79^/7) (vertex, offset from the best fit for 
clarity), 100/^/7 (edge, scaling expected for 7 > 10^), and 2OOOO/7 (face, scaling for bending of 
flat facet [ol). Inset: diametrically opposed forces applied with uniaxial symmetry. Dashed blue 
line shows 5/^/7. 



B. Spheres 



Now consider the analogous response for the case of icosahedrally symmetric triangulated 
spheres and faceted icosahedra. We consider the inverse of the susceptibility as an effective 
spring constant K = 1/x, and make the spring constant dimensionless by dividing by the 
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Young's modulus Y. We first present numerical results for symmetric forces over a wide 
range of 7 values, then in later subsections consider the limiting cases of small and large 7. 

The data was generated from a sequence of icosahedra of varying sizes and elastic prop- 
erties. We consider nonchiral icosahedra with P = Q ranging from 4 to 512 (i.e. A*"^ ranging 
from 482 to 7864322). To speed calculation, the full 120-element icosahedral symmetry 
group Yh was employed, resulting in a speedup of nearly 120 times. Beyond the buckling 
transition the five-coordinated vertices sharpen, with a radius of curvature related to the 
buckling radius 



Rb = V-fbK/Y ^ 12.4 (50) 

by a geometrical factor of order 1. In order to approximate the continuum limit we chose to 
hold kg = 1 fixed and keep kf, > 1/2, resulting in Rf, > 7.6a. 

Each structure was relaxed using a conjugate-gradient method. We found that the nec- 
essary number of relaxation steps diverges with increasing size, consistent with the 1/R^ 
vanishing relaxation rate predicted by eq. (123|) . To ensure sufficient accuracy in the suscepti- 
bility, we used 128-bit real arithmetic in the final stages of all relaxations. For P = Q = 512, 
complete relaxation requires approximately two months on a 3.0GHz Intel Xeon computer. 



For studies such as ours which seek the continuum limit, a finite element aproach |15l. Il9l. |2C 
might be more computationally efficient than our discrete mass and spring model. 

Once the structure was relaxed without applied stress, we re-relaxed with a radially 
inward force F applied symmetrically at all = 12 vertices, all = 20 faces or all = 30 
edges. The effective spring constant K = F/( was defined as the applied force F divided by 
the displacement ( of the mass to which the force was applied. We actually consider NK/Y 
because we define K as the derivative of ( with respect to all A^ simultaneous apphed forces 
F. Small applied force F = 0.001 was required to achieve linear response in cases where K 
became small. 

Figure [6] shows numerical data for symmetric forces applied to vertices, edges or faces. In 
the limit of small 7 the three data sets converge to a 7-independent value. As 7 increases, 
the vertices weaken more quickly than the faces or edges, consistent with our picture of 
the buckling transition as concentrating at the disclinations which are located at vertices. 
However, beyond 7^ the vertex stiffness falls off very slowly, while both face and edge stiffness 
continue their rapid decline. 
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1. Small 7 limit 

The following discussion first considers the limit of small 7, in which the shapes are 
nearly spherical and calculations can be done exactly. The response depends on whether the 
stress is applied in a uniaxial manner (e.g. at diametrically opposed points) or in a more 
symmetric manner (e.g. applied simultaneously at all vertices or faces or edges, or even an 
isotropic pressure). 

For an applied pressure P the deformation is purely radial, with amplitude C, as in the 
breathing mode discussed previously. This increases the energy by SnB^'^, while doing 
work AuR'^PC, against the pressure. Balancing the two yields ( = R^P/AB, susceptibility 
X = d(/dP = P? /AB and spring constant K/Y = AB/YR^. In the case of symmetrically 
applied point forces, we identify P = NF/AtxR^ yielding NK/Y = IQnB/Y = \2tx = 37.7, 
where we used eqs. (138|) and (140!) . The stiffness is independent of 7, consistent with the 
numerical result shown. 

For uniaxial stress, let the displacement at the two poles be ( and assume this displace- 
ment persists over a polar region of size d (see section 15 of Ref. |27|). The bending energy 
density ~ k{C,/(PY, and integrating over the polar region yields total bending energy 
E}, = Ki^'^/d?. Meanwhile the strain tensor Uap ^ (/R yields a total stretching energy (see 
eq. fl28l) ) Eg ~ Y{C,/ R^dP . Minimizing the sum E^ + E}, to find the optimal shape yields 

~ {k/Y)R^ and Es + E^ ^ \J kY jR. Equating this to the work done against the 
applied force, we find C ~ {Rj V kY)F and x = R/ VkY. The elastic constant K/Y ~ 1/^7) 
independent of the axis along which the force is applied, consistent with our numerical results 
(see Fig. El inset). 



2. Large 7 limit 

For 7 > 76 the radius of curvature at the icosahedron vertices quickly approaches 
Ry (eq. [50]) and remains fixed independent of the icosahedron radius R. Forces applied 
at icosahedron vertices get transfered through the curved vertex region to the fiat facets in 
a primarily longitudinal manner. According to the theory of longitudinal deformation of 



plates (see section 13 of Ref. [27|) the displacement at large distances r from the applied 



force varies as u{r) ~ (E/Y) logr/ro with vq some fixed length. Upon setting K = F/u{R) 
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and choosing tq proportional to -R^,, we find that K /Y ~ c/ log (6^/7). The numerical data 
shown in Fig. [6] fits well to this form with values c=17.3 and 6=0.79. The curve shown for 
comparison illustrates c = 15 imposing a uniform displacement for visual clarity. 

For forces applied to the icosahedron edges we expect to see the onset of ridge scaling 



behavior 
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mm, 



20, |2l|, |22| as 7 approaches 10^. Unfortunately the diverging relaxation 



time prevents us from exploring larger 7 within our current calculational method, preventing 
us from observing this behavior cleanly. We briefiy review the predictions of ridge scaling. 

Let L (which is proportional to R) be the length of an icosahedron edge. At each end of 
this edge the facets join at a fixed angle of ^ = 138.2°. At the middle the edge sags inward 
by an amount creating a saddle shaped ridge with a small radius of curvature Ri across 
the ridge and a large (and negative) radius of curvature R2 along the ridge The strain 
along the ridgeline is of order {(/L)^. Because the facets on either side of the ridge approach 
angle 9, the radius -Ri is proportional to the sag ( [l7|. Assuming that the bending and 
strain energy persist along the length L of the ridge and extend a distance Ri to either side, 
we estimate the energy as 

E = R,L[Y{C/LY + K{l/R,f] - FC (51) 

where the final term represents the action of a force F acting at mid-edge. 

Upon setting C ~ i?i and varying Ri to minimize the energy, we find, in the absence of 
force F, 



i?i ~ (/€/F)1/6l2/3 ^ v^7'^' ~ ^7"'^'- (52) 
In the presence of weak applied force F, the small radius Ri increases by an amount of order 

Ai?i ~ (53) 
V kY 

Recalling that ^ ~ i?i and converting this to an effective spring constant K = dF/ d( yields 



K/Y ~ ^JkJYL? ~ 1/v^- Indeed, the edge elasticity in Fig. [6] seems to show a crossover 
towards slope —1/2 on our log-log plot. 

Meanwhile, the icosahedron faces become almost planar in the limit of large 7. Timo- 
shenko 9j discusses the deflection of an equilateral triangular plate under load applied at 
the center. The deflection is proportional to R? j from which we conclude, using eq. ([T]), 
that KjY ~ 1/7. However, in Fig. [6] the face elasticity seems to follow a power law closer 
to -0.8 than -1. Perhaps residual stresses in the faces or on their boundaries are responsible 
for this difference. 
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V. CONCLUSIONS 

In summary, we investigated the eigenvalue spectrum of a simple mass-and-spring model 
of a virus capsid as it passes through its buckling transition. The buckling of a spherical 
shell occurs in a smooth, nonsingular fashion, in contrast to the buckling of a disclinated 
planar network. The smearing can be attributed to symmetry-breaking between the interior 
and exterior of the shell and is caused by the forbidden crossing of the buckling mode with 
a lower frequency breathing mode. 

Symmetries of the icosahedron and analogies with continuum elastic theory were used 
to classify the normal modes. Modes of full icosahedral symmetry, transforming as the 
unit representation, soften as the Foppl-von Karman number passes through the buckling 
transition. Displacements during buckling, which resemble the maturation of real virus 
capsids, can be well represented as a superposition of the two lowest icosahedrally symmetric 
modes. 

Susceptibilities to applied forces diverge at the buckling transition for planar networks. 
For spherical topology they evolve smoothly, with anomalies in the vicinity of 7;,. Suscep- 
tibility to forces applied at icosahedron vertices dominates near 76, but icosahedron edges 
and faces are much softer for large 7. In the limit of small 7 the effective spring constant 
approaches the behavior of a spherical continuum. 

Beyond the buckling transition the faces have the softest linear response, so this is where 
one might expect rupture in response to an isotropic osmotic pressure. Relative softness of 
icosahedron faces as compared to vertices has been reported experimentally in liposomes jj] . 
We verified this numerically by calculating the Qq parameter which measures the distortion 
from a sphere to an icosahedron 5|]. Below 75 isotropic pressure weakly increases the value 
of Qe, while above 7^ pressure strongly decreases Q^, bending the facets to make the shape 
more nearly spherical. 
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